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In the recent paper by Bernardini et al. [J. Comput. Phys. 232 (2013), 1-6] the discrepancy 
in the performance of finite difference and spectral models for simulations of flows with a 
preferential direction of propagation was studied. In a simplified investigation carried out 
using the viscous Burgers equation the authors attributed the poorer numerical results of 
finite difference models to a violation of Galilean invariance in the discretization and propose 
to carry out the computations in a reference frame moving with the bulk velocity of the flow. 
Here we further discuss this problem and relate it to known results on invariant discretization 
schemes. Non-invariant and invariant finite difference discretizations of Burgers equation are 
proposed and compared with the discretization using the remedy proposed by Bernardini et 
al.. 

1 Introduction 

In the recent paper [1] a possible remedy was discussed to improve the poor numerical behavior 
of finite difference simulations of turbulent flows with a preferential propagation direction. It 
was shown that the violation of Galilean invariance of the finite difference scheme is the most 
likely explanation why it is necessary to use a significantly larger number of grid points in 
finite difference calculations than in spectral methods to achieve comparably accurate numerical 
results. The recommendation given in [1] is to carry out the finite difference computations in a 
reference frame that moves with the constant stream-wise bulk velocity in the flow direction. It 
was then shown for the example of Burgers equation that the finite difference model may yield 
similar numerical results as spectral discretizations with approximately the same number of grid 
points. 

In the present paper we further discuss this problem and the remedy proposed in [1]. In 
fact, the problem found and analyzed in [1] has been investigated quite intensively in the field of 
group analysis of differential and difference equations, see e.g. [2,3,10,11,14,15,19] and references 
therein for some of the most recent results. In particular, it was established by Dorodnitsyn 
and collaborators [6, 11, 12] that it is not possible to maintain the Galilean invariance of partial 
differential equations in a finite difference model when the mesh does not move in the course 
of the numerical integration. This result qualitatively explains why the method proposed in [1] 
may work from the geometrical point of view. 

The violation of Galilean invariance of stationary discretizations can be readily shown by 
applying a Galilean boost, which in the one-dimensional case is 

(t, x, u) = (t, x + et, u + e), (1) 

where e 6 M, to the defining equation of the grid, — x™ = 0. Here and in the following, an 
upper index indicates the time level and a lower index the spatial grid point. The action of the 
Galilean transformation (1) on this grid equation yields — xf = x™ +1 — xf + e(t n+1 — t n ), 
which clearly fails to be invariant for e ^ 0. Here we assumed that all the grid points are defined 



on the same time layer, i.e. = tf = t n . It can be checked that this assumption does not 
violate the invariance of most of the equations of hydrodynamics, see also [11] for more details. 

Unfortunately, to maintain Galilean invariance it is also not sufficient to carry out the nu- 
merical simulations with a standard finite difference scheme in a constantly moving reference 
frame as proposed in [1]. It can be verified numerically that the resulting numerical solutions in 
the resting and in the convecting reference frames do not coincide, which is explicitly shown in 
Figure 1 for a FTCS discretization of Burgers equation. In this figure, we display the numerical 
solution at t = 0.5 in the resting reference frame (solid line) and in a reference frame which 
moves with constant velocity £3 = 1 (solid line with triangles) as in [1]. 




X 



Figure 1. Integration using the classical FTCS discretization of Burgers equation (2). Solid line: Original 
integration in a resting reference frame. Solid lines with triangles: Integration in a reference frame moving with 
constant velocity e = 1 as in [1]. The results in the moving reference frames were shifted back to the origin for 
proper comparison. 

Instead of using a non- invariant finite difference scheme in a convecting reference frame, it is 
therefore desirable to construct proper finite difference discretizations that preserve the invari- 
ance group of a physical differential equations. The above observation on the incompatibility 
of stationary meshes with Galilean invariance have severe consequences on the design of finite 
difference models for the equations of fluid dynamics. In fact, it renders necessary to come 
up with strategies to combine the requirement of using moving meshes (in order to preserve 
Galilean invariance) with approaches that lead to discretization schemes having good numerical 
properties, such as stability, optimal grid adaptation (e.g. equidistribution of the discretization 
error) and the possibility for a parallel implementation. From a more general point of view, it is 
necessary to bridge the fields of group analysis and numerical analysis of differential equations. 

To outline this connection for the example of Burgers equation considered in [1] is the main 
aim of the present paper. In Section 2 we discuss invariant finite difference schemes for Burg- 
ers equation. We construct three different types of invariant numerical schemes, namely La- 
grangian discretizations, invariant adaptive Eulerian schemes and invariant schemes employing 
an evolution-projection strategy. We relate these schemes to the remedy for reducing the effect 
of violation of Galilean invariance proposed in [1]. Numerical results for the different schemes 
discussed are presented in Section 3. The final Section 4 contains the conclusions of the paper. 
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2 Invariant finite difference schemes for Burgers equations 



As in [1], we introduce Burgers equation as a canonical example for high Reynolds number flows, 

u t + uu x - vu xx = 0, (2) 

where v > is the viscosity, which could be scaled to 1 by means of an equivalence transforma- 
tion. It is one of the most investigated models in the group analysis of differential equations, 
see e.g. [4,5,17]. Its maximal Lie invariance algebra q is spanned by the basis elements 

dt, d x , td x + d u , 2td t + xd x - ud u , t 2 d t + txd x + (x - tu)d u . (3) 

The associated one-parameter Lie symmetry groups are 



r 2 
r 3 
r 4 



(t, X, U) !->• (t + El, X, u), 
(t,X,U) \-¥ (t,X + E2,ll), 

(t, x, u) i->- (t, x + e 3 t, u + e 3 ), 
(t, x, u) 1 y (e 2£4 t, e £4 x, e~ £4 u), 



(4) 



T 5 : (t, x,u)\-+ ( - , - * ^ , u(l - e 5 t) + g 5 x j , 

showing that Burgers equation (2) admits time translations, space translations, Galilean boosts, 
scalings and time inversions as one-parameter symmetry transformations. 

Invariant numerical schemes for Eq. (2) have already been investigated in the literature [9- 
11,21]. The schemes constructed in these references preserve the entire five-parameter symmetry 
group G of Burgers equation. However, as was discussed in [3], it is more natural to preserve 
only those symmetries that are compatible with a particular set of initial and boundary value 
problems chosen. In the present, we focus on periodic boundary conditions. The time inversion 
T§ is then not compatible with a periodic domain as it does not map any periodic function 
u to another periodic function for £5 7^ 0. As a result we only aim to numerically preserve 
the first four symmetry transformations ri— 1^4 . These transformations form the subgroup G 1 
of the maximal Lie invariance group G of Burgers equation. It is also important to note that 
the subgroup G 1 is typical for various models of fluid mechanics. Consequently, the strategies 
discussed below are also relevant for physically more interesting higher-dimensional models of 
hydrodynamics, such as the Euler or Navier-Stokes equations. 

Arguably, the most important observation established in the field of invariant finite difference 
schemes is that it is generally not possible to maintain all symmetries of a system of differen- 
tial equations if the discretization scheme is constructed on a fixed, orthogonal discretization 
mesh [11, 15, 19]. In the case of Burgers equation, it is the presence of the Galilean transfor- 
mations T3 that prohibits the use of a fixed discretization mesh. This was explicitly shown in 
the introduction. Hence, finite difference models operating on a fixed mesh cannot be Galilean 
invariant. 

A possible remedy is to use the following expression as a discretization of Eq. (2) 

< +1 ~ < , (n _ -A <K ~ <! _ ^ _ = ,5, 

v ' x i+l x i-l \ x i+l x i x i x i-l / 

where xf = (x™ +1 — xf)/ At. Applying the transformations T1-F4 one readily verifies the invari- 
ance of this discretization. The reason for this scheme being invariant is that the introduced 
grid velocity xf transforms as xf — > xf + £4 under the action of the Galilean transformation and 
the additional term involving £4 is exactly compensated by the Galilean transformation of uf. 
In other words, introducing a moving mesh into the discretization of Burgers equation restores 
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the Galilean invariance of the finite difference scheme. See [11, 14, 18, 19] for further details on 
the systematic construction of invariant finite difference discretization schemes. 

In [3] it was shown that the above discretization can be interpreted as a discretization of (2) 
in terms of computational coordinates, i.e. in a coordinate system that remains fixed in the 
presence of grid adaptation. To accomplish this transformation, one sets x = x(6,£), where 
9 = t and £ = £(t,x) is the spatial computational coordinate. Transforming (2) to the (#,£) 
coordinates and discretizing the result using a FCTS scheme leads to (5). 

The question remaining is how to determine the grid velocity xf, which involves the yet un- 
known location of the grid points on the subsequent time layer n+1. As pointed out in [19], there 
are two main strategies to find x™ : The first is to construct a grid equation that is invariant 
under the same symmetry group as the discretization of the physical differential equation. The 
second method is to regard the grid adaptation as unconstrained from the symmetry require- 
ments imposed by the physical differential equation, i.e. to use a non-invariant grid equation. 
We will mostly focus on the first method here as it is geometrically more grounded. 

In the first method we require the grid equation to be invariant under the same symmetry 
group as is the discretization of the physical differential equation. In the present case, this 
amounts to constructing a grid equation that is invariant under the transformations T\— One 
simple possibility is to take 



as this equation is obviously Galilean invariant and also does not violate the remaining trans- 
formations from G . This choice boils down to discretizing Burgers equation in Lagrangian 
coordinates, i.e. the grid velocity equals the physical velocity. Indeed, Lagrangian discretization 
schemes are among the earliest examples of invariant discretizations for hydrodynamical equa- 
tions, see e.g. [11]. The problem with Lagrangian discretizations is that one in general does 
not have proper control over the evolution of the grid points. This can be a severe problem, 
especially in the multi-dimensional case, where grid points can concentrate in certain regions, 
deteriorating the local resolution of the scheme in areas away from these concentration regions. 

Perhaps numerically more satisfying are grid equations that couple the evolution of the grid 
to the development of pronounced features in the numerical solution, i.e. to use proper grid 
adaptation strategies. Linking grid adaptation to the construction of invariant discretization 
schemes proved relevant in the numerical investigation of blow-up problems, see e.g. [7,8,13]. A 
possible way to realize an invariant adaptive grid is based on the equidistribution principle for 
a monitor function p, 



which plays a central role in the construction of r-adaptive numerical schemes in one space 
dimension. See again [8, 13] and references therein for an extensive discussion of the concept 
of equidistributing meshes. An invariant equidistributing mesh is obtained by discretizing (7) 
in a G 1 -invariant way. The general feasibility of this approach depends on the structure of the 
symmetry group one aims to preserve [3] but in general it can be realized for the symmetry 
groups one usually encounters in hydrodynamics. The basis for this approach is to choose a 
proper monitor function p, that leads to a form of (7) that is invariant under the symmetry 
subgroup G 1 of Burgers equation and then to discretize this expression in an invariant way. A 
(^-invariant monitor function is 



which coincides with the arc- length function for a = 1. The reason for including a generic a in 
this expression is that the term u x is not scale invariant, i.e. it transforms as u x = e~ 2£3 u x and 
thus, as it stands, p is not scale invariant. Extending the scalings of (t, x, u) to an equivalence 



xf - u? = 



(6) 



{pxi)i = 0, 



(7) 



p = \A + au l 
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transformation involving a by adopting the transformation rule a = e 2e3 a then indeed leads to 
a G 1 -invariant function p. It should be stressed though that no such extension of the Galilean 
transformation to an equivalence transformation is needed to guarantee the Galilean invariance 
of the resulting form of the equidistribution principle (7). 

Using the modified arc-length weight function, a possible (^-invariant discretization of (7) is 



n /-I , U i+l 



(8) 



which can be solved using a relaxation scheme, such as e.g. GauB-Seidel iteration to obtain 
and hence to complete the invariant numerical scheme (5). More sophisticated ways to 
solve (7) are conceivable as well and could be used to improve the quality of the resulting 
adaptive discretization scheme. 

A further possibility for the construction of invariant numerical schemes is to invoke an 
evolution-projection strategy. This idea was put forward in [16, 20] for the non-invariant dis- 
cretization of advection equations and extended in [2] to find an invariant discretization of the 
linear heat equation. The main approach in the evolution-projection strategy is to use the in- 
variant numerical scheme and the invariant mesh equation only for a single integration step and 
to use a projection operator (i.e. an interpolation) to map the solution from the off-grid points 
back to the initial, uniformly spaced mesh. If the interpolation is done in an invariant way, i.e. 
the interpolation used preserves the invariance (sub)group of the system of differential equations 
being discretized, then the entire discretization procedure becomes invariant. The advantage of 
this approach is that moving meshes can be completely avoided. 

In the present case of Burgers equation (2) we observe that classical interpolation schemes 
such as linear, quadratic or cubic spline interpolation already preserve the invariance subgroup 
G 1 . This means that we can use the aforementioned interpolations to re-map the solution -u™ + 
defined at the points back to x™ +l S { x f} 5 without breaking the invariance of the scheme. 

We show this explicitly for linear interpolation here, which is defined as 

_ ..n+1 , u i+l a i (-n+1 _ „n+l\ _ n+1 n+1 n+1. ~n+U 

u l x i ) — u i ' n+1 _ n+1 \ X i X i > ~ M x i > x i+l ' u i ' u i+l ' x i I 

X i+1 X i 

for the interpolation of values lying within the interval [x™ +1 , x"^ 1 ]. Then, for transforma- 
tions of the form r^— IT^, we obtain that 

?(% a+1 \ — r(z- n+l ^■ n + 1 1 -; n + 1 l. n+1 \ — r,,t^ n + l \ _ r(r r n +^ ™«+! „,"+! 

u \ x i I L -\ X i i x i+\i u i ,U i+1 ,X i ) — U{X { ) L.\X i , X i+1 , U i ), 

which proves the invariance of linear interpolation under transformations from G . Similarly, 
invariance of quadratic and cubic spline interpolation can be shown. 

Regarding the use of non-invariant grid equations, which is the second possibility to complete 
the description of the scheme (5), in principle all choices excluding = xf are admissible. 
We will now return to the remedy proposed in [1]. The proposed approach falls into the category 
of non-invariant grid equations. In that paper, it was suggested to use a reference frame moving 
with the (constant) bulk velocity of the flow. In the case of Burgers equation, the authors set 

= x n % + cAt, (9) 

where c = const. It is readily verified that this grid equation is not Galilean invariant, as 

~n+l _~n_ = x n+l _ % n _ ^ _ 

If c is the bulk velocity of a flow, one may indeed expect that the numerical results obtained 
from scheme (5) with grid equation (9) are better than for scheme (5) on a stationary grid. 
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More precisely, Galilean invariance could be restored by extending the transformation T4 to c by 
setting c = c + £4. This extension is justified in case c is related to u, which is the main reason 
why the remedy proposed in [1] may work. In other words, reference frames moving with a 
constant velocity could be made Galilean invariant in the sense that Galilean transformations 
have an extension to equivalence transformations for such reference frames. A similar extension 
of the transformation T3 to c is necessary to incorporate the scale invariance. 

For the sake of convenience, we summarized the characteristics of the different schemes dis- 
cussed in the present section in Table 1. 



Table 1. Different numerical schemes for the viscous Burgers equation (2). 



Numerical scheme 


Grid equation 


Galilean invariance 


Grid spacing 


Finite differences 


None 


Not invariant 


constant 


Lagrangian 


Lagrangian grid movement 


fully invariant 


variable 


Eulerian adaptive 


Equidistribution principle 


fully invariant 


variable 


[1] Bernardi et al. 


Constant grid movement 


semi-invariant 


constant 


Evolution-proj ection 


Lagrangian grid movement 


fully invariant 


constant 



3 Numerical results 

In this section we present some numerical results obtained from the invariant discretization 
schemes introduced in the previous section. For all the experiments, we use u(0, x) = sin(:r) 
as the initial condition on a 27r-periodic domain and the integration is carried out up to time 
t = 0.5 for v = 0.1. Unless otherwise stated, we use N = 64 grid points and fix the time step 
with At oc /i 2 , where h is the mean grid spacing over the domain. 

In Figure 2 we carry out numerical integrations for three invariant numerical schemes for 
Burgers equation based on (5) and employing different grid equations: (i) The fully Lagrangian 
scheme uses the Lagrangian grid equation (6), (ii) in the invariant Eulerian adaptive scheme 
the grid points is determined from the invariant discretization of the equidistribution 

principle (7) and (iii) in the evolution-projection scheme we use the Lagrangian grid equation (6) 
and use quadratic interpolation to re-map the off-grid points to their original location at the 
previous time step. 

For all schemes we numerically verify Galilean invariance, i.e. each of the pairs of integrations 
in a resting and a constantly moving coordinate system yields visually the same numerical 
solution for (i) the Lagrangian scheme (dashed-dotted line, £3 = and dashed line, £3 = 1), 
(ii) the Eulerian adaptive scheme (solid line with diamonds, £3 = and solid line with crosses 
£3 = 1) and (iii) the evolution-projection scheme (solid line with squares, £3 = and solid line 
with circles £3 = 1). Moreover, it is seen from Figure 2 that all three schemes yield approximately 
the same numerical solution. We aim to detail and analyze this result further now. 

The difference between the three invariant numerical schemes for Burgers equation is the 
invoked grid equation. For the sake of a clearer presentation, we depict the grid spacing Axf = 
^f+i — as a function of the location of the grid points x^ at the final integration time t = 0.5 
in Figure 3. For the classical, non-invariant finite difference scheme the spacing is by definition 
constant (solid line). For the Lagrangian discretization (dashed line) the location of the grid 
points depends on the solution itself and therefore one does not have control over the local 
resolution. It is a mere consequence of equating the physical velocity and the grid velocity. In 
the adaptive Eulerian scheme (dashed line with diamonds) we observe a proper concentration of 
the grid points along the building shock, which follows from using the equidistribution principle 
and the arc-length monitor function. Away from the steepening front, the points remain quasi- 
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Figure 2. Integration of Burgers equation (2) using the three invariant discretization schemes introduced in 
Section 2. Dashed-dotted line: Lagrangian scheme in the resting reference frame. Dashed line: Lagrangian 
scheme in a reference frame moving with constant velocity £3 = 1. Solid line with diamonds: Eulerian adaptive 
scheme in the resting reference frame. Solid line with crosses: Eulerian adaptive scheme in a reference frame 
moving with constant velocity £3 = 1. Solid line with squares: Evolution-projection scheme in the resting 
reference frame. Solid line with circles: Evolution-projection scheme in a reference frame moving with constant 
velocity £3 = 1. 



— Classical finite differences 

- - ■ Lagrangian scheme 

0— Eulerian adaptive scheme 
a— Evolution-projection scheme 




Number of grid point 



Figure 3. Grid spacing A:r™ = x" +1 — x™ at final time t = 0.5 for the four discretizations of Burgers equation 
shown in Figures 1 and 2. 



equally distributed. By construction, the grid points in the evolution-projection method (solid 
line with squares) are again equally spaced. 

To estimate the overall accuracy of the different schemes, in Figure 4 we display the pointwise 
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Figure 4. Pointwise difference of the numerical and the exact solution at final time t — 0.5 for the four discretiza- 
tions of Burgers equation shown in Figures 1 and 2. 



differences of the numerical solutions and the exact solution of the chosen initial value problem 
it(0, x) = sin(x) for Burgers equation, which is 



u e (t,x) 



2v 



Y.7=i a jJ e vtj2 sinjx 



■\oo 

J _i 



i 



2tt 



1 



a = — I e-V-^'Mdx, a j>0 = - I , 
2^ Jq vr J 



2n 



-(1-cobx)/(2v) 



cos jx dx. 



It can be seen from Figure 4 that all the numerical schemes achieve a comparable accuracy. The 
overall -errors of the schemes for the runs depicted in Figure 4 are summarized in Table 2. 

Table 2. Zoo-errors for the various discretization schemes for the viscous Burgers equation (2) with N = 64 grid 
points. 





classical FD 


Lagrangian 


Eulerian adaptive 


Evolution-projection 


\\E\\i 

II II too 


2.53 • 1(T 3 


1.69 • 1(T 3 


2.50- 1(T 3 


2.63 • 1(T 3 



The convergence rates of the different schemes in the Zoo-norm are depicted in Figure 5 using 
N G {4, 8, 16, 32, 64, 128, 256, 512} grid points. It is seen from this plot that all schemes yield 
comparable errors. The overall convergence rates demonstrate that all three types of invariant 
schemes for Burgers equation introduced in this paper are asymptotically of second order and 
that the different invariant discretization strategies do not alter the accuracy of the underlying 
approximations. More details on this can be found in [2]. 



4 Conclusion 

In the present paper we have revisited the problem recently pointed out in [1] that classical 
finite difference discretizations of the governing equations of hydrodynamics violate Galilean 
invariance. We have discussed three possible ways of constructing Galilean invariant finite 
difference schemes for Burgers equation: (i) Lagrangian discretization schemes, (ii) Eulerian 




Figure 5. Convergence plots for the four numerical schemes presented above for TV 6 {4, 8, 16, 32, 64, 128, 256, 512} 
grid points. Solid line: Classical finite difference scheme. Dashed-dotted line: Lagrangian scheme. Solid line with 
diamonds: Eulerian adaptive scheme. Solid line with squares: Evolution-projection scheme. 



adaptive discretizations and (iii) discretizations using an evolution-projection strategy. These 
three approaches can be readily adapted to the two- and three-dimensional Euler or Navier- 
Stokes equations. The approach proposed in [1] leads to a semi-invariant scheme in that the 
authors use the discretization on a moving coordinate system but employ a non-invariant grid 
equation by moving the grid points with a constant velocity (e.g. the bulk velocity) rather than 
with the actual physical velocity. If properly done, this approach can indeed reduce the effects 
of the violation of Galilean invariance in classical finite difference discretizations, while still not 
a fully Galilean invariant schemes. 

All of the invariant discretization approaches presented above have their advantages and 
disadvantages. Purely Lagrangian discretization schemes are not in widespread use as such 
schemes usually lead to a strong concentration of grid points, leaving other regions of the 
domain poorly resolved. Moreover, in multi-dimensional cases of interest in hydrodynamics, 
Lagrangian schemes can lead to tangled meshes. At the same time, Lagrangian schemes are 
good in that they are able to preserve sharp interfaces within a fluid. The Eulerian adaptive 
approach is attractive because it provides a natural way to link the problem of finding invariant 
discretization schemes to the properties of the numerical solution at each time step. At the 
same time, the computational overhead required to efficiently generate the meshes can be a 
crucial factor determining the feasibility of the adaptation methodology, especially for multi- 
dimensional systems. Finally, in situations where adaptive numerical schemes are not desirable, 
the evolution-projection strategy based on the Lagrangian grid equation (or any other invariant 
grid equation one is able to find) and a simple (but invariant) interpolation is a possible way to 
maintain Galilean invariance in a finite difference scheme while still being able to operate on a 
fixed, uniformly spaced mesh. The drawback of the evolution-projection approach is that it also 
requires an additional operation, the interpolation, which may cause additional computational 
overhead compared to, for instance, the Lagrangian scheme. 

Applied to Burgers equation, the three invariant discretization methodologies yielded numer- 
ical schemes that are asymptotically second order accurate. Compared to the standard FTCS 
scheme, which is also second order accurate, one thus gets as a bonus the preservation of a geo- 
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metric feature of Burgers equation, namely a subgroup of its symmetry group. The preservation 
of this symmetry subgroup can be a crucial factor for problems for which it is inevitable to carry 
out the simulations in a moving reference frame. The construction of higher order invariant 
numerical schemes is the subject of current research and will be reported in the future. 
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